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The phase diagram of a simple model with two patches of type A and ten patches of type B 
^-1— ( (2yllO-B) on the face centred cubic lattice has been calculated by simulations and theory. Assuming 

that there is no interaction between the B patches the behavior of the system can be described in 
terms of the ratio of the AB and AA interactions, r . Our results show that, similarly to what happens 
for related off-lattice and two-dimensional lattice models, the liquid-vapor phase equilibria exhibits 
reentrant behavior for some values of the interaction parameters. However, for the model studied 
here the liquid- vapor phase equilibria occurs for values of r lower than | , a threshold value which 
was previously thought to be universal for 2AnB models. In addition, the theory predicts that below 
■| (and above a new condensation threshold which is < |) the reentrant liquid- vapor equilibria 
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is so extreme that it exhibits a closed loop with a lower critical point, a very unusual behavior in 
^ \ single-component systems. An order-disorder transition is also observed at higher densities than 

the liquid-vapor equilibria, which shows that the liquid-vapor reentrancy occurs in an equilibrium 
region of the phase diagram. These findings may have implications in the understanding of the 
condensation of dipolar hard spheres given the analogy between that system and the 2AnB models 
^ . considered here. 
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I. INTRODUCTION 



Advances in the fabrication of nanometer-to-micrometer sized particles enable tailoring their size, shape and intcr- 
" actions, but their organization into complex structures remains a challenge. Self- assembly is an appealing route of 
this bottom-up approach as the structure of the clusters is tunable through the anisotropy of the particle shapes and 
• • ■ interactions. In addition, the strongly anisotropic interactions prevent the clustering that drives condensation and 
. £h i have been shown to lead to novel macroscopic behavior, including empty liquids, optimal networks and equilibrium 
^ . gels±r£. 

$_i ' Patchy particle models with dissimilar patches (A and B) were introduced in this context^ and allow a unique 
5^ ■ control of the effective valence through the temperature T. In three-dimensional (3D) off-lattice models consisting of 
particles with two types of patches, A and B, with no interaction between the B patches, the topology of the liquid- 
vapor diagram is determined by the ratio between the AB and the AA interactions, r = sab/^aa- As r decreases 
in the range | < r < i, the low-temperature liquid- vapor coexistence region also decreases^. The binodal exhibits 
a reentrant shape with the coexisting liquid density vanishing as the temperature approaches zero&£. Below r = | 
condensation is no longer observed, and above r = j there is no reentrant behavior—. 

Both the scaling of the vanishing critical parameters and the reentrant phase behavior are predicted correctly by 
Wcrtheim's thermodynamic first-order perturbation thcory£i2r— . The theory also reveals that the reentrant phase 
behavior is driven by the balance of two entropic contributions: the higher entropy of the junctions and the lower 
entropy of the chains in the (percolated) liquid phase, as suggested a decade ago on the basis of a hierarchical theory 
of network fluids^. 

In a previous paper we considered the 2A2B model consisting of particles with four patches, two of type A and 
two of type B, on the square lattice and investigated its global phase behavior by simulations and theoryi^. We have 
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set the BB interaction to zero and calculated the phase diagram, as a function of r. We found that, in the same 
range of parameters as in the 3D off-lattice models, the liquid-vapor diagram exhibits a reentrant shaped, with a 
region where the system exhibits empty fluid behavior—. In addition, below r = ^ condensation ceases to exist, while 
the reentrant regime disappears for r > 5, in line with the results for 3D off-lattice models and the predictions of 
Wcrtheim's theory^. When r < | the gain in entropy resulting from the AB bonds does not balance the loss in 
energy of the favored A A bonds, in line with the simulation rcsults^^. This led us to suggest that the thresholds 
predicted by Wcrtheim's theory are exact and universal, i.e. independent of the dimensionality of the system and of 
the lattice structure^. 

Computer simulations of the phase diagram of these models become more and more demanding as r decreases. One 
of the reasons is the rapid increase in the size of the voids of the empty coexisting liquid, as the temperature decreases. 
Simulations of larger and larger systems are thus required to obtain reliable results, rendering the computation of 
phase equilibria in the neighboorhood of r = | prohibitive^. 

In this paper we present results that clarify the role of the | threshold and the universality of the empty fluid 
regime reported earlier. We consider a model consisting of particles with twelve bonding sites ("patches"), two of 
type A and ten of type B, on the face centered cubic lattice, and investigate its global phase behavior by simulation 
and theory As before we set the interaction between the B patches to zero. The potential energy of this 2AnB class 
of models, is given by: 

U = -caaNaa - zabNab = - (Naa + tNab) £aa; (1) 

where eab and cab are positive, and Naa and Nab are the number of AA and AB bonds, respectively. 

The model is a 3D counterpart of the 2A2B model on the square lattice^. We develop efficient sampling algorithms 
to simulate the empty liquid at low temperatures and observed liquid- vapor coexistence in systems with r < i. 

We establish a new non-universal threshold r m < i for liquid-vapor equilibrium (LVE) and Wertheim's first-order 
perturbation theory suggests the existence of a new regime where the reentrant behavior is so extreme that the low 
temperature binodal closes at a lower critical point. Using an asymptotic expansion of the theory, we calculate the 
new threshold and show that the closed loop regime exhibits some degree of universality as it depends only on the 
number of B patches. The conditions to observe the closed miscibility loop are met in lattice models with a large 
number of B patches, such as the 2A\QB lattice model, but are unlikely to occur in continuum models with the 
same number of patches^. The existence of a threshold r m < i for LVE and its degree of universality have physical 
relevance, as it has been used to address the liquid-vapor condensation of dipolar hard spheres^. 

In more recent work, the absence of liquid-vapor coexistence of dipolar fluids was related to ring formation^, a 
feature which is not described by Wertheim's first-order perturbation theory on which the reported thresholds are 
based. A related study of the generic phase diagram of 2A4B models on the triangular lattice with r > |, reports that 
orientational correlations between the A patches that promote ring formation have a profound effect on the phase 
equilibria^. No phase coexistence is observed when short rings are formed. Closed miscibility loops are found if 
larger rings are formed while the usual reentrant behavior is observed if no rings are formed. Somewhat surprisingly 
the same regimes are reported in this work based on Wertheim's first order perturbation theory, which does not 
account for ring formation, where r is the control parameter. The orientation of the AA bonds that promotes ring 
formation, on the triangular lattice, has an effect on the topology of the phase diagram which is similar to the effect 
of decreasing the AB interaction in a system without rings but with a large volume available for the formation of AB 
bonds. Beyond this observation, the relation between the generic phase diagrams of 2AnB systems, with and without 
rings, is an important open question that will be addressed in future work. 

The remainder of this paper is arranged as follows: In Sec. [II] we describe the model. In Sec. Mil we describe the 
Monte Carlo techniques used to compute the phase diagrams. In Sec. IIVI we present the simulation results while 
in Sec. [V] we carry out the theoretical analysis. Finally, in Sec. I VII we discuss these and previous results and the 
perspectives for future work. 

II. THREE DIMENSIONAL MODEL 

We consider a face-centered-cubic (FCC) lattice. Sites on the lattice can be either empty or occupied by, at most, 
one particle. The particles carry twelve patches: two of them of type A, and ten of type B. The patches on each 
particle are oriented in the directions linking the site with its twelve nearest neighbors (NN). The angle between the 
two A patches is 180 degrees and the line between these patches defines the orientation of the particle. On the FCC 
lattice, the particles are oriented along one of the q = 6 equivalent directions Si = s(f\) = 1,2, ••• ,q and thus a 
lattice site has q + 1 possible states (9 = 6 orientations plus the additional empty state s» = 0). The grand canonical 
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Hamiltonian can be written as H = U — fj,N, with /i being the chemical potential and N the number of occupied sites: 

M 

N = "£[l-5 0tSi ); (2) 

i=l 

where 8 is Kronecker's delta function and M the total number of lattice sites. The potential energy IA can be written 
as: 

M q 

U = ^2^2V2[s(f f i ),s(r i + a k )] (3) 

i=l fc=l 

where V2 is the interaction between pairs of NN sites, which depends on the states of the sites and the direction of 
dk (the vector linking the two sites). The pair interaction is defined as: 
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In short, interactions are defined between two NN sites, i and j, if both are occupied, with a magnitude that 
depends on the types of patches on each site pointing to the other one. We take e = caa as the energy scale and in 
line with previous work set: €bb = 0. 

The 2A1QB model with ess = and r = eAB/t < \ may exhibit two phase transitions (as in 2Di^) which are 
interpreted as a liquid- vapor transition, and an order disorder transition, the latter occurring at a higher density when 
both transitions occur at a given temperature. 



III. SIMULATION METHODS 



We have adapted the methodology used in 2Dii to the 3D lattice model. In addition we have built up cluster 
algorithms to enhance the sampling procedures. In 3D the order-disorder transition was found to be discontinuous, 
which simplifies the calculation of the phase diagram. We used cubic cells of different sizes and periodic boundary 
conditions. The number of sites of a given system is M = 4L 3 , with L an integer. 



A. Liquid-vapor transition 



We have considered different values of r in the range [0.30,0.50]. The basic protocol to compute the liquid- vapor 
equilibria (LVE) was as follows: First we used Wang-Landau multicanonical (WLMC) procedures^—, supplemented 
with a finite-size scaling analysi o 19 i 21 ~ — to compute the value of the chemical potential of the transition at some 
subcritical temperature, To. This temperature is chosen not too close to the critical temperature of the LVE to 
find a clear separation between the two phases, but not too far to avoid the sampling difficulties that appear at low 
temperatures. For different system sizes, L, we computed the value of the chemical potential that maximized the 
density fluctuations of the system, /x e (L,To); this is considered the finite-size estimate for the LVE at temperature 
To- Then, the value in the thermodynamic limit, fi e (To), is obtained by fitting the results to: 

ti e (L,T ) = fi e (T ) + aL~ d ; (5) 

where d is the spatial dimensionality of the system, d = 3. The range of L required to get a reliable estimate of ^t e (To) 
depends dramatically on the value of r. As one reduces r, larger values of L are required. The results were checked 
by performing a fully independent estimate by means of thermodynamic integration techniques using relatively large 
system sizes. Details of the implementation of these techniques will be described later in the paper. 

Once we have estimated a reference point for the LVE: [To,/i e (To)] a Gibbs-Duhem integration procedure was 
carried out to draw the binodal lines, using a methodology similar to that described in previous paper o 14 i 24 ' 25 . The 
only relevant difference is that in the grand canonical simulation of each phase, we incorporate collective moves via 
the cluster algorithms described in the Appendix. 
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1. Wang-Landau multicanonical methodology 



The basic strategy of the multicanonical procedures is to sample in a single run the properties of the system with 
different numbers of particles (N). The procedure can be related to a simulation in the grand canonical ensemble 
(GCE). In both cases the probability of a given configuration, S^r (with N occupied sites), can be written as: 

P(S N ) = lu q (N) exp \-pU (Sjv)] , (6) 

where f3 = (fc^T) -1 , T being the temperature, and fee the Boltzmann's constant. In the GCE simulation the weighting 
factor is taken to be: u>o(N) oc cxp (JjfiN) whereas in the WLMC method, ujq(N) is chosen to obtain a flat histogram of 
the density P(N) ~ 1/(1 + N max - N min ); WN e [N min7 N max ], which in practice implies u) (N) cx exp [f3F(N, M, T)], 
F being the Hclmholtz thermodynamic potential. In order to compute ujq {N) in the WLMC procedure an equilibration 
scheme based on the Wang-Landau strateg y 26 i 27 is carried out, where the values of u>o(N) can vary through the first 
part of the simulation^. Once wo(iV) is obtained, the equilibrium simulation (fixed wo(iV)) is run, and values of the 
different properties are computed as a function of N. These results can be used to determine phase equilibria, and 
using appropriate reweighting techniques to locate the critical points. 

The WLMC simulations include two types of moves: insertion and deletion of particles on the lattice. In most cases 
(r > 0.31) we used a simple non-biased algorithm: particles to be removed, and the position and orientation of the 
inserted particles were chosen at random (among non-occupied lattice sites). Taking into account detailed balance^, 
the acceptance criteria for these moves can be written as: 

A(S w+1 | Sjv ) = mm|l,— _ e j\T+T J & 



^(S,- 1 |S,)= m i n {l ) ^e--- I ^ T ^}; (8) 

where AJ7i ns and At/del are respectively the variations of the energy when inserting and deleting one particle. The 
factor q in Eqs. ([7][H]) takes into account the q possible orientations of one particle, whereas the M — N like terms are 
the number of empty sites where the insertions can be attempted. 

For the systems with the lowest values of r (r < 0.31) (and at low temperatures), the previous scheme was found 
to loose efficiency. Then, in addition to the random insertions and deletions described above, we include biased 
insertion/deletion moves. In the biased moves only insertions (deletions) that lead to the formation (destruction) 
of, at least, one AA bond are considered. The insertions are exclusively tried in empty sites to which at least one 
non-bonded A patch is pointing. We choose at random with equal probabilities one of those non-bonded A patches, 
and perform the trial insertion of the new particle at the site to which the chosen patch is pointing to, with the 
orientation that guarantees the formation of an AA bond with that patch. Deletions are carried out following the 
symmetric move, i.e., we choose at random, with equal probability, one of the A-patches that participates in an AA 
bond, and the particle to which the A patch points to is removed. Taking into account the conditions of detailed 
balance, the acceptance criteria of these moves can be written as: 

A ; (f5jv+i Sjv) - mm <M, — e — = — — (9) 

I uj (N) 2Naa{Sn+i) J 



A^ (S N ^\S N ) = min (l, ^ e -^ |W | (10) 

I u (N) Nao{Sn-i) J 

where Maa is the number of AA bonds in the system, and Mao is the number of A patches that point to an empty 
site. Before attempting a MC step we choose with equal probabilities whether to use fully random or biased moves. 



2. Critical points 



The critical points of the LVE were estimated using the Wang-Landau multicanonical algorithm. After preliminary 
short calculations to locate approximately the critical temperature we run, as above, simulations for different system 
sizes. By storing histograms of the mean values of the potential energy: < U > and its square < U 2 > as a function 
of the number of particles TV we can apply a reweighting scheme to the results to estimate the thermodynamics at 
temperatures close to that where the multicanonical simulation is carried outJ^ 
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We proceed to calculate the pseudo-critical points (pi L , 7c ), with pi L ^ = p e (Tc L \L). We compute the cumulants 
gi(L,T, p e ) defined by (74 = rr^/m 2 ,, where are the k th order moments of the density distribution function 
P(p\T,fx e ,L): 

m k =< (p - pf >; (11) 
and p is the average of the density. The pseudocritical point satisfie a 14 ' 19 i 23 ' 29 : 

g 4 (L,T<; L \^)=gi c) ; (12) 

(c) 

where g\ is a constant that depends on the universality class and on the boundary conditions. The liquid- vapor critical 

point of 2AnB models is expected to be in the 3D Ising universality class, with g± ~ 1.604M. The pseudocritical 
densities are taken as: 

p c (L)=p(L,T( L \fxi% (13) 



The estimates for the critical properties in the thermodynamic limit, (T c , p c ), are then obtained by extrapolation of 
the pseudocritical values using the scaling law o 23 i 29 . 

T C (L) - T c oc (14) 
p c {L) - p c oc L' 1 '". (15) 



where and v are critical exponents 
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3. Gibbs-Duhem Integration 

In order to compute the binodals of the LVE we used a version of Gibbs-Duhem integration (GDI)2i in the 
GCE 14 ' 24 i 25 . GDI requires the knowledge of an initial point (T , p ) on the phase coexistence line. Then, the binodals 
are obtained by solving the differential equation: 

dfl = (E - d T. (16) 

p \T TAN J v ' 

The equation is solved numerically using finite intervals, and a fourth-order Runge-Kutta procedure. The A's represent 
the difference between the mean values of the properties (U or N) in the two phases at equilibria, computed by 
simulations of both phases for the same system size. The sampling of the simulations was enhanced by incorporating 
cluster moves as described in the Appendix. We used L = 32 as the system size for the integrations, except at low 
temperatures where L was increased to L = 64. 



B. Order-disorder transition 



The order-disorder transitions were computed using thermodynamic integration techniques which give the value of 
the chemical potential at coexistence, no, at a given temperature. This can then be used as the initial coexistence 
point to calculate the binodals using Gibbs-Duhem integration. 

The initial coexistence point is obtained by choosing a value of the chemical potential sufficiently high to ensure that 
the system is fully occupied at T = 0, and computing the equation of state p(po, T) for two sequences of temperatures, 
one starting at infinite temperature (disordered phase), and the other starting at T ~ (ordered phase). In these 
limits the grand potential per site is known exactly. The coexistence temperature is determined by the equality of 
the grand potential of the two phases. 



1. The ground state 

In the ground state (GS), the stable configurations minimize the grand canonical Hamiltonian: % = U — pN . The 
minimum energy of the model with N occupied sites occurs when all the A-patches are bonded through AA bonds, 
i.e. the potential energy is Uqs = —Ne, which implies: 



U GS (N) = (-e-p)N; 



(17) 
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It follows that the equilibrium state corresponds to the empty lattice when /i < — e, and to the full lattice with M 
AA bonds when \x > — e. The fully occupied ground state is highly degenerate^ as a large number of configurations 
are compatible with Uqs — —Me. In general, the order of a fully occupied configuration may be described through 
the order of a set of planes [i,j, k]: Particles in one plane have the same in-plane orientation but different planes may 
exhibit different (in-plane) orientations. There are two types of such planes [1, 1, 1], and [1, 0, 0]. In [1, 1, 1] planes the 
sites form a triangular lattice, one site has six NN on the plane, and there are three in-plane orientations. In [1, 0, 0] 
planes the sites form a square lattice, one site has four NN on the plane and there are two in-plane orientations. 
Taking into account the boundary conditions, the degeneracy of the GS at full occupancy is then: 

O ^ (3 x 4 L ) + (4 x 3 L ) (18) 

It follows that the entropy of the ground state scales linearly with L (S oc L) and that the entropy per site vanishes 
in the thermodynamic limit. 

The grand canonical partition function S at very low temperatures and at full occupancy (fx > — e) is: 

So(M, fi) = O cxp [+/3Me + /3Mfx) , (19) 

from where the grand potential can be easily calculated using ^ 

$ = -fc B TlnS, (20) 

It follows that at low temperatures in the thermodynamic limit: 

4>=^- = -e-LL- (T^0:iJ>-e:L^w). (21) 



2. The high temperature limit 
In the limit of high temperatures, /3e a( g — > 0, the grand canonical partition function can be written as: 

s = E ( n ) w 

N=0 ^ ' 

where we took j3fx finite. The factor q N accounts for the q possible particle orientations. Using elementary combina- 
torics we find: 

E=(l + qe^) M , (23) 

which, as /3 — >• at finite fi, yields 

LnS = Mm(l + g). (24) 
Then, the grand canonical potential per site at high temperatures is given by: 

^ = -ln(l + g); (T^oo); (25) 



3. Computing the initial order-disorder coexistence point 

Given the grand canonical potential of one point in each phase, the order-disorder transition is located using 
thermodynamic integration. In the Grand Canonical ensemble the variation of the thermodynamic potential at 
constant volume (M) is given by: 

UN H 
d (/?<£) = — d/3 - — d (/3a0 = — d/3 - pTidp, (26) 

where U is the potential energy, defined by Eqs. ([3]|4]), and 77 = N/M. This may be used to calculate the thermody- 
namic potential 4>{T) of the ordered phase as a function of temperature by thermodynamic integration at constant 



d(^) = h(T)d(3, 



(27) 
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FIG. 1: Grand canonical potential of the ordered and disordered phases for r—0.32 at fj,/e = —0.9. The order-disorder phase 
transition occurs at the crossing of the two branches. 



where we defined h = %/M. In order to calculate the grand canonical potential of the ordered phase Eq. [37] has to 
be rewritten to avoid divergences at low temperatures. First we add and subtract the value of the integrand at zero 
temperature and at full occupancy (which is ho = — e — fio): 

d{/3<j>{T))= hod/3 +[h(T)- h ] d/3, (28) 

and then change variable from /? to T. The grand canonical potential of the ordered phase is given by: 

MT) = Ph h{T ^; h ° dT>, (29) 

where the integrand is well behaved over the domain of integration, approaching zero as T — > 0. 

Similarly, the grand canonical potential <p(T) of the disordered phase may be computed by integrating Eq. 1271 fusing 
(3 as the integration variable) from the infinite temperature limit at constant volume and fj, = jio'- 

= - ln(l + q) + / h{/3'W, (30) 
Jo 

The order-disorder transition is discontinuous and for sufficiently large systems exhibits hysteresis. This simplifies 
the calculation of the transition temperature at fixed [i as it is possible to evaluate H(T) for each phase in a range 
of temperatures around coexistence. The two branches of $(T) (ordered and disordered) cross at the transition 
temperature (see Fig. [1]). We computed the transition temperatures at \iojt = —0.90 for several values of r: 0.00, 
0.32, and 0.40. In all cases we found negligible size effects and used L = 24 , and L = 32 as typical system sizes. 

Having obtained the transition temperature at fx = jio, GDI integration is used to calculate the coexistence lines at 
different temperatures, as was done for the liquid- vapor equilibria. At low temperatures we used large system sizes 
L = 64, otherwise L = 32. 



C. Consistency checks 

The calculations of the liquid- vapor and order-disorder transitions were checked by performing simulations with two 
fully independent programs written by two of the coauthors of this paper. In one of the MC codes cluster sampling 
techniques were included whereas in the the other (control program) only simple single site moves were included. The 
control code was used to calculate the liquid-vapor and order-disorder transition using thermodynamic integration. 
These calculations validate both the enhanced sampling techniques and the proper coding of the Wang-Landau and 
cluster moves in the GCMC/GDI simulations. 

The calculation of the order-disorder transition using the control code was carried out using the thermodynamic 
integration technique described in the previous section, whereas the calculation of the liquid-vapor equilibria by 
thermodynamic integration is described briefly in what follows. 
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For the vapor phase, the grand canonical potential was obtained integrating along an isotherm from very low values 
of the chemical potential (or densities), where the system behaves as an ideal gas, to the chemical potential of interest 
using: 

P<t>{n) = m^deai) -P I r)(»'W (31) 

where fXideal represents a value of the chemical potential low enough so that the behavior of the fluid can be considered 
ideal. The grand canonical potential of the ideal gas can be easily calculated using: 

M(n) = -p P (32) 

and the ideal gas equation to obtain: 

/34>(fJ.zdeal) = -rj (33) 

The free energy of the liquid was calculated through an integration path starting at the high temperature limit, 
where the free energy is calculated as in the previous section (Eq l25p . The integration was performed in two steps: 
first, we integrate from infinite temperature to the temperature of interest keeping the chemical potential fj, constant 
(Eq. [27]) and, second, we integrate along an isotherm from /j, to the chemical potential of interest (Eq. [3Tjl . As usual 
we ensure that there arc no phase transitions along the chosen thermodynamic path. 

Analogously to the calculation of the order-disorder transition the liquid-vapor coexistence is obtained by calcu- 
lating, at the given temperature, the chemical potential where the grand canonical potentials of the two phases are 
equal. 

All the simulations performed with the control Monte Carlo code considered systems with L =24. The results 
obtained using the two different codes and methodologies are found to be consistent. 



IV. RESULTS 



In Figure [5] we plot the results for the liquid- vapor equilibria, including the estimates for the critical point (these 
are also given in Table U). The general trend is as expected from earlier wor k 6 ' 7 ' 14 : As r decreases the critical point is 
shifted to lower densities and lower temperatures. In addition for r < 1/2 the LVE is reentrant, with densities of the 
liquid phase at coexistence decreasing on cooling at low temperatures. The LVE binodals were computed using GDI 
with system size L = 32 (i.e. M = 131 072) for r > 0.30 (except at the lowest temperatures, where L = 64 (M = 
1 048 576) was used to avoid interconversion between the two phases). For r = 0.30 larger systems, L = 64, were 
required to obtain reliable results at all temperatures. In all cases, the GDI fails at sufficiently low temperatures, due 
to the rapid growth of the typical size of the voids in the (emptying) coexisting liquid as reported previously^. As in 
other 2AnB models 6 -^ ' 14 ' 18 the binodal at a given r encloses the binodals at smaller values of r. 

The most remarkable result, however, is the clear evidence of LVE for systems with r < |. As mentioned in the 
Introduction, earlier theoretical predictions based on Werthcim's first-order perturbation theory set the threshold for 
liquid-vapor coexistence at r = ^ in line with simulation results for 2AnB models 6 -^. The finding is also relevant in 
a wider context, as this threshold was used recently to address the liquid- vapor condensation in other systems that 
form branched chains at low temperatures, in particular dipolar hard spheres^ 6 -. We will return to this discussion 
later, after the theoretical analysis described in Sec. [V] 

Lattice models allow the precise location of not only the liquid- vapor transition, but also the order-disorder transition 
that occurs at higher densities (the lattice analogue of the fluid-solid transition of off-lattice models) . As in the 2D 
lattice, we find that the order-disorder transition occurs always at a higher density (higher chemical potential), in the 
temperature range accessible to simulations (see Fig. [3]^. 

Interesting scaling features are revealed by plotting the LVE and the order-disorder binodals for different values 
of r as functions of the scaled temperature, t = fcgT/[(l — 2r)e] (see Fig. [3]). First, in the limit of full occupancy 
the order-disorder transitions collapse into a single point. This is an exact result, and the observed collapse is just a 
check of the consistency and accuracy of our simulation protocols. In addition, the order-disorder transition exhibits 
a very weak dependence on r; i.e. the lines for the order-disorder transition for r = 0.32, and r = 0.40 are almost 
indistinguishable. By contrast, the order-disorder transition of the SARR model (r = 0), deviates clearly from the 
previous ones. Finally, as in the T-77 representation, the binodal for a given r encloses the binodals for lower values 
of r. 
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11 

FIG. 2: Liquid- vapor binodals, for different r = exs/e. Large open symbols mark the critical points. Continuous lines represent 
GDI results, the points on these lines mark the portions that were computed using larger system sizes (L = 64). Dashed lines 
are included to connect the GDI results with the critical point estimates. 



r 


T* 

± c 






0.300 


0.1163(9) 


0.105(4) 


-1.0361(10) 


0.305 


0.1277(3) 


0.1284(10) 


-1.0522(4) 


0.310 


0.1371(4) 


0.148(2) 


-1.0679(6) 


0.320 


0.1518(2) 


0.175(2) 


-1.0961(3) 


0.330 


0.1644(1) 


0.1993(10) 


-1.1245(2) 


0.350 


0.1860(1) 


0.2331(2) 


-1.1808(2) 


0.400 


0.2302(1) 


0.2848(3) 


-1.3228(1) 


0.450 


0.2689(1) 


0.3136(3) 


-1.4684(1) 


0.500 


0.3050(1) 


0.3310(2) 


-1.6169(1) 



TABLE I: Estimates of the critical points for different values of r 



V. WERTHEIM'S THEORY FOR 2AnB LATTICE MODELS 



The free energy per particle, within Wertheim's first order perturbation theory, for a homogeneous system of 
particles with 2 patches of type A and n patches of type B is^£, 

Pf = Pfref + 2(lnX A - + n(\nX B - ^), (34) 

where f re f is the free energy per particle of the reference system and X a the fraction of unbonded patches of type 
a. The laws of mass action that relate X a , the density rj and the temperature T arc (when there are no BB 
interactions)^, 

X A + 2 V A AA X 2 A +m 1 A AB X A X B = l, (35) 
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FIG. 3: Binodals for the liquid-vapor and order-disorder transitions, for different r = cab /e, as functions of a rescaled 
temperature. Filled symbols mark the portions of the curves that were computed using larger system sizes (L=64). 



X B + 2 V A AB X A X B = 1. (36) 
As for the 2A2B model on the square lattice, the reference system is an ideal lattice gasi^, and thus, 

[3f ref = In r] + ln(l - n). (37) 

V 

The A a/ 3 are integrals of the Mayer functions of two patches a and j3 on two different particles, over their positions 
and orientations, weighted by the pair distribution function of the reference system. In continuous systems the A Q /j 
are calculated from, 

A a/ 3 = J drxdr 2 j ' (^i(^29 re f(ri,f 2 )f a p(f 1 ,f2,oji,uJ2), (38) 

where a is a particular patch on particle 1 and /3 is a particular patch on particle 2, r, refers to the position of particle 
i and ui to the orientation of the particular patch on particle i that is being considered; the factor l/(47r) 2 takes 
into account that there is no preferred orientation for the position of the patches on the particles surfaces. Finally, 
fafsiri, Tii wi, LU2) is the Mayer function of the interaction potential between patches a and (3, and g r ef(ri, is the 
pair correlation function of the reference system^. 

The calculation of A Q( g on a lattice with coordination number z and particles with z patches, (as in the 2AnB 
models considered here and ir*i£) is carried out by discretizing Eq. [3^1 



^ AT M z z 

P = TT^Y1 H ^2 fap(ii,h,h,h)gref(ii,i2)- (39) 



Mz 

ii=l 12 = 1 h = l (2=1 



Here the- patch a is on particle 1 and the patch fi on particle 2; ii and i 2 represent the positions of particles 1 and 2, 
respectively; M is the total number of lattice sites (equivalent to the volume); the factor l/z 2 accounts for the fact 
that there is no preferred orientation for the patches, as z is also the number of different orientations of a given patch. 
The integers l\ and / 2 run over the possible orientations of each patch. For the potential described in Sec. jTl] the 
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Mayer function f a p(ii,i2,li,l2) is non zero when: (a) 1 and 2 are NN and (b) l\ and I2 are such that patches a and 
(3 are properly oriented along the interparticle direction. Using g re f = 1, Eq. [35] is simplified to, 

A Q(3 = v a/3 [exp(/3e Q(3 ) - 1] , (40) 

where v a p, the volume of a bond between patches of type a and /3, is 

v a!3 = 1/z. (41) 

For the 2AnB models under consideration, the bonding volume v a p is independent of the types of patches and is 
related to the number of B patches z = n + 2. 

The liquid-vapor binodals calculated using Wertheim's theory, for models with n = 10, i.e. the models simulated 
in Sec. IIV1 are plotted in Figure 21 for several r. As expected&ii^, Wertheim's theory predicts, in agreement with 
the simulation results, the reentrance of the liquid binodal for r < |. The phase diagrams were calculated for models 
with r > 0.305: below this value no liquid-vapor coexistence was found. 

To clarify this surprising behavior, in view of previous results^— we solved the set of equations, (^j^Jz^ = and 

(frp^) = ^> wn i cn determine the critical density and temperature as a function of the parameters of the lattice 2AnB 

model, namely r and n (or the coordination number z). The results for 71 = 2, n = 4, n = 6 and n = 10 (corresponding 
to square^, triangular or simple cubic, body centered cubic, and face centered cubic lattices, respectively) arc plotted 
in Figs. [5] and[6l 

The results reveal that Wertheim's theory predicts a phase behavior that is strongly dependent on the values of n 
and r. For n < 4 we recover the threshold reported earlier: no critical point for r < | and one critical point otherwise. 
For larger n, however, three distinct regimes are possible depending on the value of r: one critical point for r > -|, 
no critical point for r less than a threshold r m (< |) and two critical points for the range of r between this threshold 
and 3. 

A deeper understanding of these results, which contrast with the simpler picture reported earlier— is obtained 
through the asymptotic expansion of the free energy Eq. [2H in the limit of strong AA bonding (i.e Xa ~ 0) and low 
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FIG. 5: Critical temperature as a function of r = 6ab /e, for several values of n. Full lines: critical temperature calculated 
from Wertheim's theory Eq. 1341 dashed lines: critical temperature calculated from the asymptotic expansion of Wertheim's 
theory, Eqs. ((43j) and (|44|) . 



densities^ — . Using this expansion we find for the pressure p, 

r/2 nA AB 



V2A 



AA 



V2A 



B2V 2 , 



AA 



(42) 



where B2 is the second virial coefficient of the reference system (B^ = 1/2 for the ideal lattice gas). Using Eq. [42] the 
critical temperature T c is found to satisfy, 



G{eAB,T c 



exp 



- 1 



exp 



= C, 



(43) 



where C = ^vabY a ' ^ or ^AnB models under consideration, the constant C may be evaluated using Eq. l4Hand 
B 2 = 1/2, 



C 



2{n + 2) 2 



(44) 



The asymptotic critical density rj c is then given by, 



Vr 



nA AB , c 



(45) 



where Aab,c is the value of Aab at T = T c . In Figures [5] and [6] the asymptotic critical temperature and density are 
plotted as functions of r for several values of n. As expected, the asymptotic results describe those of the full theory 
at low temperatures, but qualitative agreement is obtained at intermediate temperatures. Thus, the phase behavior 
of the 2AnB model can be described by analyzing Eq. 031 
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FIG. 6: Critical density as a function of r for several values of n. Full lines: critical density calculated from Wertheim's theory 
Eq. 1341 dashed lines: critical temperature calculated from the asymptotic expansion of Wertheim's theory, Eqs. (|43[) , (|44|) and 
l|35). 



For r > i the function G(cab, T c ) is a monotonic decreasing function of T c , with limits lim^-^o = 00 and liniT c -->oo = 
0. However, for r < g, this function has a maximum which decreases with decreasing r; is limited by < G(eAB, T c ) < 
1 and vanishes in the limits, lim-j^o = and lirriT c ->oo = 0. The analysis of the critical behavior of 1AnB models 
is done most simply by considering the cases C > 1 and C < 1 (which correspond, Eql4"3] to n < 4 and n > 4, 
respectively, for integer values of n). 

C > 1: If r > |, Eq. [43] has a unique solution, and therefore, there is a single critical point with critical temperature, 

(46) 



InC 

If r < I then G(e J 4B,7c) < 1 and Eq. [43] has no solution. Therefore, there is no critical point. 
This is the case considered in previous worker— , as the models analyzed previously have C > 1. 

C < 1: If r > I, Eq. [43] has one solution and there is a single critical point. If r < |, we define r m as the value of r 
where the maximum of G, G max , is equal to C: G max {r mi T c ) = C Then, two cases have to be distinguished: 
If r, n < r < i, G max (eAB, T c ) > C, and Eq. [43] has two solutions; therefore, there are two critical points; If 
r < r m , Eq. SSJhas no solutions, and there is no critical point. 

In Figure [7] we plot a phase diagram with two critical points, for models with n = 10, and r = 0.31 and 0.32. 

The liquid vapor coexistence of these models is between a low density, high energy and high entropy phase, formed 
by short chains, and a high density, low energy and low entropy phase (network liquid) formed by long chains connected 
by AB bonds, or junctions^. It has been shown^ that this coexistence is only possible when a decrease in the entropy 
of chains (or AA bonds) upon condensation, is balanced by an increase in the entropy associated with the junctions 
(or AB bonds) . For systems with C > 1 and r < | , the increase in the entropy of the junctions is no longer sufficient 
to balance the loss in the entropy of the chains^. 

Systems with C < 1 have not been discussed earlier as the continuum^ and the lattice^ models investigated 
previously belong to the class C > 1. The 2A10B lattice model investigated in this paper belongs to the class 
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FIG. 7: Phase diagrams with two critical points, for n — 10 and two values of r < 1/3, calculated using Wertheim's theory Eq. 
1341 The circles represent the location of the critical points. These binodals are the same as those of figure [4] for the indicated 
values of r = tAs/e, but are represented in a logarithmic scale in the density and extended to lower temperatures, to put 
highlight the lower critical points. 



C < 1 and thus exhibits different critical behavior. In these models the balance of entropies may occur at values 
of r < |. This may be rationalized by considering the physical meaning of C. The entropy of one bond is the 
logarithm of the volume available (on one particle) to form that bond^ 4 -. Then, InC = \n( 2 ^ A ) — 3 ln( B ) is the 
difference between the entropy of one ^4^4 bond and (three times) the entropy of one AB bond. For the 2AnB model 
InC = ln(^g) — 31n(^p5), and as n increases so does the entropy of the AB bonds. Therefore, when C < 1 the 
entropy of junction formation increases, in such way that it can balance the decrease of entropy of the chains, for 
values of r < 3 . 



VI. DISCUSSION OF THE RESULTS 



Despite the challenges posed by the simulations of the phase diagram of empty fluids at low temperatures, the 
results for the 2A10B lattice model confirm the features of the LVE reported for 3D off-lattice 6 -^ and 2D lattice 
2AnB^ models. The variation of the critical densities and temperatures with r follow the expected behavio r 6,7 ' 14 . In 
addition, we computed the order-disorder transition that occurs at higher densities, confirming that the low density 
liquid phase is thermodynamically stable as in the 2D model^. 

We have, however, found an unexpected result: LVE for systems with r < ^, by contrast to previous simulation 
results and the theoretical analyses based on Wertheim's theory 4 ^ 7 - as well as an earlier prediction based on a hierar- 
chical theory of network fluids 1 ^. The threshold r — ^ results from an asymptotic expansion of Wertheim's first-order 
perturbation theory, which assumes that the constant In C is positive as described in Sec. |Vl This is certainly the 
case for 2AnB models on and off-lattice if the number of B patches is not too large. For lattice models, however, 
the bonding volume compatible with a single bond per patch assumed by Wertheim's theory is much larger than in 
similar off-lattice models (the exclusion of the second particle being guaranteed by the lattice structure) and thus In C 
can become negative. In this case, Wertheim's first-order perturbation theory and its asymptotic expansion in the 
limit of strong AA bonding predicts indeed the possibility of LVE for r < |. The theoretical analysis also predicts 
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that in this regime the reentrancy of the liquid-vapor binodal is extreme in the sense that the system exhibits a low 
temperature critical point. The theoretical prediction is then that when InC is negative (large values of uvab) 2AnB 
models exhibit a closed miscibility loop in a range of r < |. There is also a new threshold, which depends on the 
number of B patches, below which the closed miscibility loop vanishes and where there is no condensation. 

Previous simulation results on and off lattice were compatible with the original ^ threshold, and in line with the 
theoretical results for positive In C&£^. 

The closed miscibility loop predicted for systems with negative InC has not been confirmed by simulations, since 
the density and temperatures at which they occur are beyond the current simulation techniques. 

In related work, a 2AAB 2D lattice model with \ > r > \ was shown to exhibit the usual reentrant behavior when 
the position of the A patches prevents the formation of rings, a closed miscibility loop when the orientation of the A 
patches promotes relatively large rings and no phase coexistence when the orientation of the A patches promotes short 
rings-^. The topology of the phase diagram of this 2B4A lattice model with r > ^ changes as the orientation of the 
A patches changes (promoting the formation of rings) in a fashion that resembles the behavior of the 2A10B model 
as r decreases. Although the physics may be related a detailed, quantitative and qualitative, analysis is required in 
order to investigate the analogies in the driving mechanisms of the different transitions. 

Along these lines, recent work for a model with A patches addressed quantitavely the competition between ring and 
chain formation, within an extension of Wertheim's first-order perturbation theory and by simulation^. An extension 
of this approach to 2AnB models and the calculation of the corresponding phase diagrams is a challenging task that 
will be addressed in future work. Likewise new simulation algorithms will be developed to confirm the presence of 
closed miscibility loops, in systems with no rings, as predicted by Wertheim's first-order perturbation theory as well 
as the new thresholds, r < |, for models with negative InC. 

The degree of universality of the new thresholds is also an important open question, in general, and in the context 
of the condensation of dipolar hard-spheres. 
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VII. APPENDIX: CLUSTER ALGORITHMS 

In this appendix the cluster algorithms that we used in the Grand Canonical ensemble simulations of the GDI 
method are described. The two cluster moves defined here do not include, in general, the whole set of sites of the 
lattice, but only those sites with two of the possible values of Sj. In practical terms we can classify the cluster moves in 
two types: Moves that change the number of occupied sites: Cluster N-sampling, and moves in which the orientation of 
some of the sites can change: Cluster Orientation-sampling. A full derivation of the procedures might be cumbersome, 
so we will just include the steps and considerations required to understand the recipe of the algorithms. 

A. Cluster N-Sampling 

In these moves we consider only empty sites and occupied sites with one (chosen at random) of the six possible 
orientations, s r . These sites are named active sites. We classify as passive (or blocked) those sites k with Sk ^ s r and 
Sfc 7^ 0. Passive sites are not modified in these moves, and play the role of an external field. Taking into account the 
values of the interaction parameters, and in particular that ebb = 0, one occupied active site interacts with another 
occupied active site if (and only if) both sites are NN in the s r direction. Using this definition of active sites the 
system can be equivalently seen as a collection of one dimensional rows of sites (with PBC), i.e. ID lattice gas models 
under the influence of external fields. The terms of the Grand Canonical Hamiltonian that deppend on the active 
sites can be written as: 

A A 

H A = -e PiPj-YsP^V + trN^ + ermii)} (47) 

<ij>s r i 
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where the first sum on the right hand side is carried out exclusively over pairs of active sites which are NN along the 
direction s r . The second sum includes only active sites. The variables pi take the values for empty sites and 1 for 
occupied sites. N\(i) is the number of A patches of the site i that points to a NN passive site, and A^i) is the number 
of A patches belonging to a NN passive site of i that point to site i. Through the change of variables pi = (1 + Oi)/2; 
"Ha, may be written as an Ising-like Hamiltonian: 

A A A 

H A -H = - e - £ a i a j -^Y,^-Y,^i N ^) + {^- e i) N ^ ( 48 ) 

<ij>s T i i 

where Ho includes the terms that do not depend on the state of the active sites. The new variables <x; can take the 
values ±1. On the right hand site of Eq. (|48p the first term is the Ising-like interaction, the second plays the role 
of a global external field, and the last one includes the local external fields that depend on the configuration of the 
passive sites. 

Within this representation of the interactions of the active sites, it is straighforward to build up a cluster algorithm 
following the Swcndsen-Wang procedure 3 -^ and its extensions in the presence of external fields^. The recipe of such 
an algorithm goes as follows: (1) Generate bonds between pairs, of active sites which are NN in the s r direction, 

and that fulfill Sj = Sj with probability: 

Ba = 1 - exp [-e/(2k B T)\ ; (49) 

(2) Consider separately each one of the clusters of active sites defined by the previous bonds. Taking into account the 
effect of the external fields given in Eq. (|4"B1 , the new configuration is generated by assigning to all the lattice sites in 
the cluster either s = (i.e. a = —1); or s = s r (i.e. a = 1) with probabilities A(c,s) (where c is the index for the 
c-th cluster) fulfilling: 



A(c,s r ) (p + e (r- l/2)eA/i(c) + erAf 2 (c) 
AM) = 6XP 1 n * (c) + fcsT 

where rij(c) is the number of lattice sites in the cluster c, A/i.(c) is the number of A patches in the cluster c that point 
to a NN blocked site, and M%(c) is the number of A patches lying at blocked sites that point to a NN site in the 
cluster c. 




B. Cluster orientation sampling 



In these cluster moves two of the possible orientations, s a , s&, of the particles are chosen as active directions, whereas 
the remaining four directions and the empty sites are classified as passive (or blocked) directions. In the moves only 
active sites can modify their states (from s a to sj, and vice versa). Therefore the number of occupied sites will remain 
constant. Notice that the interaction between two active sites that are NN through a passive direction is equal to €bb 
independently of their respective orientations. In addition, the interaction between two NN sites, one being active 
and the other passive can be modified in these moves only if they are NN through an active direction. 

From these features of the interaction potential, it follows that the only relevant interactions in the proposed 
restricted sampling are those that take place between NN sites through the active directions. As a consequence the 
system can be treated as a set of independent layers with the topology of the square lattice (defined by the two active 
orientations), that can contain active sites and two types of passive sites: empty and blocked sites. The relevant terms 
of the potential energy on each layer take the form: 

H A = -e S Si)3j S Suai:j - cr^25 Styazk (1 - 6 Sk:0 ) ; (51) 

<ij> [ik] 

where the <5's represent Kronecker delta functions, < ij > indicates pairs of NN active sites on the square lattice; 
[ik] stands for pairs of NN with i and k being respectively an active and a passive site; and ai m is the index of the 
direction fi m . 

Now, we describe the strategy to generate the cluster algorithm. In previous paper s 14 ' 38 we showed how the lattice 
patchy models defined on the square lattice at full occupancy can be mapped onto the two-dimensional lattice gas 
model. It can be shown that it is also possible to carry out a mapping when some of the sites are blocked, the main 
difference being that the effect of the blocked sites enters as a local external field^ This mapping can be obtained 
through the plaquette procedures used in previous papers. Taking into account plaquettes of four site s 24 ' 38 , and 
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defining on each plaquette Ha as the sum of the potential energy contributions involving at least one active particle, 
it can be shown that Ka can be written in terms of a Potts-like interaction as: 

p v P 

h A = h - K ^2 S si, Sj + K\ ^2 ^si,a ik (1 - <$ Sfc ,o) + Kio ^2 S s ^ aik S Skfi ; (52) 

<ij> [ik] [ik] 

where the superscript p over the sums indicates that only interactions between sites belonging to the plaquette are 
considered, ho is just an additive constant which deppends on the configuration of the passive sites in the plaquette, 
but not on the state of the active sites, and finally K, Kiq and Ko deppend only on the energy parameters of the 
patchy model: e and r. Since Ha can be written as one half of the sum of the plaquettes energies Ha, we find 

H a -H = -KY, S Si<Sj + K^Ss,,^ (1 - 5 Skfi ) + K 10 £ S Si , aik 6 Sk , (53) 

<ij> [ik] [ik] 

where Hq is an additive constant which does not deppend on the configuration of the active sites, K = (1/2 — r)e, 
K\ = K, and K w = e/2. On the right and side of Eq. ([5"3"]l : the first term is a Potts q = 2 interaction; the second 
term includes the effective interaction of active sites with their NN occupied passive sites (on the square lattice), 
whereas the last term represents the effective interactions of active sites with their NN passive empty sites (on the 
square lattice). The last two terms can be seen, as before, as local external fields. 

Once the interaction between active sites has been described in terms of the Potts model, it is straighforward to use 
the same strategy described for the cluster N-sampling. Following Ref. (H3), the algorithm recipe is: Pairs of active 
sites, being NN (on the chosen square lattice) that fulfill Sj = Sj are bonded with probability: 

By = 1 - cxp l-K/(k B T)} . (54) 

These bonds define clusters of active sites; and the orientation of each of the clusters in the new configuration is 
chosen from the active directions: s a = s a , with probabilities: 



A(c, s a ) oc exp 



-JV (0) (-«)^-* W ("«)^^ 
Zkb-L kb-L 



(55) 



where N^°\s a ) is the number of A patches belonging to the cluster that point to empty blocked NN sites when the 
particles of the cluster are oriented in direction s a , and N^(s a ) is the number of A patches belonging to the cluster 
pointing to occupied blocked NN sites when the cluster is oriented in direction s a . 
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